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tension  specimen  geometry.  Two  different  crack  lengths  were  considered.  HY- 
100  steel  parameters,  previously  characterized  for  large  strain  deformation, 
were  utilized.  Additional  material  responses,  based  on  the  HY-100  nonlinear 
response  but  with  different  yield  strengths  and  ductilities,  were  also 
considered  to  parametrically  assess  material  effects  on  crack  growth,  A 
debonding  algorithm  was  employed  to  produce  crack  growth  by  nodal  relyase 
when  local  material  conditions,  satisfying  a  specified  local  fracture  criterion, 
were  met.  Material  fracture  and  crack  growth  were  treated  as  dependent 
variables  of  the  analysis,  generating  crack  growth  in  discrete  increments. 

>2-'  S'; 

The  results  of  this  parametric  computational  study  demonstrated  crack  growth 
of  up  to  sixty-seven  percent  of  the  initial  crack  length  over  the  total  number 
of  load  increments  allowed  in  each  case.  The  current  crack  tip  energy  density 
histories  exhibited  piecewise  smooth  behavior,  consistent  with  the  changing 
crack  tip  changed  position  during  crack  growth.  For  the  materials,  geometry, 
loading  and  finite  element  meshes  considered,  the  anticipated  load  versus 
displacement  discontinuities,  corresponding  to  the  discrete  crack  growth 
increments,  were  approximately  five  percent  or  less  of  the  current  load.  The 
relative  loads  and  crack  growth  sustained  by  each  specimen  were  seen  to 
depend  on  the  elastic  stress  and  strain  of  the  material  and  the  local  fracture 
toughness  of  the  material. 


OD  Form  1473,  JUN  SC  (Reverse)  security  classification  of  this  page 


ii 


CONTENTS 


l.  INTRODUCTION  .  1 

n.  REVIEW  OF  PAST  INVESTIGATIONS  .  1 

m.  OBJECTIVES  .  7 

IV.  DEBONDING  SIMULATIONS  .  7 

V.  RESULTS  .  25 

VI.  SUMMARY  .  41 

REFERENCES  .  44 


COMPUTATIONAL  INVESTIGATION  OF  LOCAL  MATERIAL 
STRENGTH  AND  TOUGHNESS  ON  CRACK  GROWTH 

I  Introduction 

Stable  crack  growth  in  ductile  materials  is  a  fundamentally  challenging 
problem.  Material,  geometry  and  load  combine  to  determine  fracture  tolerance 
or  fracture  criticality  of  a  structure.  Significant  crack  growth,  if  defined  as 
crack  growth  between  minimum  detectable  crack  size  and  maximum  allowable 
crack  size  prior  to  crack  instability,  may  be  a  large  fraction  or  even  many 
multiples  of  the  minimum  detectable  crack  size.  Structures  made  from  modem 
ductile  alloys  are,  in  principle,  capable  of  sustaining  such  significant  crack 
growth.  Prediction  of  significant  crack  growth  is  necessary  to  determine  the 
rate  at  which  crack  instability  is  approached  and  if  function  will  be  impaired 
prior  to  crack  instability.  Accurate  theoretical  and  computational  predictive 
methodologies  are  therefore  required  to  develop  a  fundamental  understanding 
of  the  crack  growth  process  and  to  utilize  ductile  materials  effectively  for  this 
class  of  materials. 

In  this  investigation,  a  typical  range  of  local  fracture  toughness  values 
are  parametrically  examined  in  the  context  of  nonlinear  finite  element 
analysis  of  stable  crack  growth.  A  debonding  algorithm  is  employed  during 
the  analysis  to  produce  crack  growth.  The  debonding  algorithm  incorporates  a 
specified  local  fracture  criterion  which  implements  a  nodal  release  scheme  to 
produce  crack  growth.  Material  fracture  is  treated  as  a  dependent  variable  of 
the  analysis,  generating  crack  growth  in  discrete  increments.  These 
increments  satisfy  local  fracture  and  equilibrium  conditions  as  imposed  by  the 
material,  geometry  and  load.  Of  particular  interest  were  local  material 
responses  at  the  crack  tip,  global  load  versus  crack  mouth  opening 
displacement  and  sensitivity  of  crack  growth  to  material  strength  and  fracture 
toughness. 

II  Review  of  Past  Investigations 

A  wide  range  of  investigations  have  been  performed  to  examine  stable 
crack  growth  in  ductile  materials.  Finite  element  analysis  continues  to  be  the 
major  computational  technique  used  for  quasistatic  ductile  material 
simulation.  This  is  particularly  true  for  materials  sustaining  large  plastic 
deformation  requiring  incremental  solution.  The  strong  influence  of  material 
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plasticity  on  the  stable  crack  growth  process  has  been  seen  in  studies  of  crack 
behavior  in  ductile  materials.  Table  1  summarizes  finite  element  investigations 
which  have  addressed  stable  crack  growth  problems.  Common  issues 
encountered  during  the  modeling  process  have  included  (i)  the  selection  of  a 
fracture  criterion  for  use  or  study,  (ii)  the  material  model  used,  (iii)  the 

material  selected,  (iv)  the  cracked  specimen  geometry  and  (v)  the  numerical 
crack  growth  scheme.  Each  of  these  issues  is  briefly  reviewed.  Specific 
references  [1-21]  for  each  issue  may  be  found  in  Table  1. 

Fracture  criteria  have  been  postulated  on  the  basis  of  local,  crack  tip 
and  global  phenomenon.  Global  approaches  were  developed  as  extensions  to 
the  energy  balance  and  stability  models  of  elastic  fracture.  Criteria  based  on 
crack  tip  parameters  treat  the  asymptotic  behavior  of  stress  and  strain  field 

solutions  in  the  immediate  vicinity  of  the  crack  tip  as  the  dominant  factor  in 

crack  growth.  Analytically  path  independent  integrals  also  related  the 
asymptotic  field  amplitude  to  far  field  quantities  amenable  to  common 
experimental  measurements.  Local  criteria  have  been  developed  to  explicitly 
address  the  various  responses  a  fixed  material  composition  and  processing 
history  can  exhibit  under  different  local  loading  conditions. 

As  fracture  criteria  address  smaller  scale  phenomenon,  the  analytical, 
computational  and  experimental  requirements  for  implementation  increase. 
The  continued  and  growing  interest  in  local  criteria  is  due,  in  part,  to  the  wide 
range  of  inelastic  material  responses  observed  and  the  cumulative  influence 
such  responses  have  on  quasistatic  crack  growth.  The  ability  to  address  the 

necessary  requirements  to  implement  such  criteria,  often  theoretically  and 
computationally  formidable,  are  gradually  being  met. 

A  number  of  incremental  plasticity  models  are  currently  in  use. 
Incremental  plasticity  models  incorporating  coupled  damage  variables  have 
also  recieved  considerable  attention  in  recent  years.  (See,  for  example,  the 
reviews  by  Chaboche  [22],  Murakami  [23],  Nquyen  [24]  and  Sarfarazi  [25].) 
Materials  selected  for  -investigation  generally  reflect  a  compromise  based  on 
relevance,  laboratory  and  testing  capabilities,  computational  capabilities  and 
technical  objectives.  Consideration  of  alloys  exhibiting  extensive  ductility, 
such  as  effective  logarithmic  strains  of  over  one  hundred 
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percent,  has  become  possible  and  increasingly  common.  Developments  in  both 
computational  hardware  and  software  have  facilitated  these  investigations.  For 
both  types  of  models,  regardless  of  the  specific  formulation  chosen  to  describe 
the  inelastic  behavior,  accurate  assessment  of  material  parameters  is  required 
for  the  full  range  of  material  response  leading  to  fracture. 

Three  major  approaches  have  been  taken  to  model  crack  growth  using 
finite  element  techniques.  These  are  element  removal,  remeshing  and  node 
release.  Each  has  been  used  with  some  success  for  a  variety  of  different 
problems.  A  brief  discussion  of  these  techniques  is  included  here. 

Element  removal,  which  associates  fracture  with  element  stiffness  loss, 
is  relatively  simple  to  implement.  The  effective  loss  of  material  at  the  crack  tip 
significantly  alters  the  crack  tip  geometry.  For  materials  failing  in  a  diffuse 
manner,  such  as  concrete,  element  removal  has  been  successfully  used.  For 
metals,  however,  the  crack  tip  remains  relatively  distinct.  Severe  deformation 
and  material  separation  remain  focused  at  the  crack  tip  throughout  crack 
growth. 

Remeshing  intermittently  recreates  the  computational  model  and 
allocates  the  necessary  resources,  in  terms  of  element  size  and  type,  to  the 
current  critical  regions  of  the  model.  For  linear  elastic  fracture  problems,  the 
process  is  straightforward.  For  fracture  problems  involving  incremental 

plasticity,  mesh  redefinition  accounting  for  prior  material  history  produces 

significant  computational  overhead  for  the  simulation  as  remeshing  occurs. 
Large  gradients  also  complicate  remeshing. 

Node  release,  which  generates  new  surface  in  a  finite  element  model, 
remains  a  popular  technique  to  model  fracture.  Two  element  edges,  initially 

constrained  to  identical  displacements,  are  allowed  to  separate  by  releasing  the 
constraints  and  nodal  forces  which  hold  the  elements  together.  Crack  surface 

is  formed  in  a  discrete,  element  by  element,  succession  of  node  releases.  If 
crack  growth  is  considered  as  an  independent  variable  for  purposes  of 
analysis,  the  nodal  release  is  provided  as  a  specified  nonlinear  constraint.  If 
crr.ck  growth  is  not  specified  a  priori  but  is  to  be  predicted,  as  is  physically  the 
case,  a  criterion  for  crack  growth  must  be  implemented.  The  criterion  assesses. 
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on  an  ek.^cnt  by  clement  basis,  whether  specified  conditions  for  crack 
growth  are  satisfied  throughout  the  analysis.  If  conditions  for  crack  growth 
are  indicated  then  node  release  is  invoked.  In  practice,  a  crack  path  through 
the  finite  element  mesh  is  specified  to  facilitate  implementation  and 
computational  efficiency. 


Ill  Objectives 

In  this  investigation,  the  influence  of  local  material  behavior  on  the 
fracture  of  ductile  materials  is  addressed.  A  simple  local  fracture  criterion  is 
implemented,  as  determined  from  experimental  measurements  and 
computational  simulations  of  material  test  specimens,  relevant  to  stable  crack 
growth.  Crack  growth,  generated  as  a  result  of  local  material  response,  was 
computationally  simulated  using  nonlinear  finite  element  analysis 
incorporating  a  debonding  algorithm. 

With  this  modeling  strategy,  a  parametric  study  of  debonding  sensitivity 
to  material  and  geometry  was  conducted.  Effects  of  material  yield  strength  and 
local  fracture  toughness  were  assessed.  Different  crack  lengths  were 
considered  to  examine  relative  crack  growth  rates  predicted  for  different 
geometries. 


IV  Debonding  Simulations 

A  series  of  debonding  simulations  were  undertaken  to  investigate  the 
two-dimensional  Mode  I  crack  growth  of  HY-100  steel  for  two  different  crack- 
length-to-specimen-width  ratios  (ao/w)  of  0.136  and  0.250  in  a  modified 
compact  tension  specimen  geometry.  Additional  material  responses  based  on 
the  HY-100  data,  but  with  different  yield  strengths  and  ductilities,  were  also 
considered.  The  two  dimensional  simulations  provided  quantitative  insight  into 
both  the  physical  and  computational  aspects  of  stable  crack  growth  in  ductile 
alloys.  The  modeling  techniques  also  serve  as  logical  prerequisites  to  full  three 
dimensional  modeling  in  the  future.  Specific  details  of  the  simulations  are 
discussed  below. 
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Specimen  Geometries  and  Loading  -  Modified  compact  tension  specimen 
geometries  (Figure  1)  were  selected  after  a  series  of  preliminary  geometry  and 
finite  element  mesh  refinement  studies.  Applied  stress  loading  at  the  upper 
left  region  of  the  specimen  was  used  to  generate  crack  opening  deformation. 
This  loading  produced  stable  crack  growth  and  was  a  computationally  efficient 
method  of  approximating  the  effect  of  standard  compact  tension  specimen 
loading  holes.  (Models  with  loading  holes  in  a  mesh  of  this  type  were  found  to 
induce  plastic  hinge  formation  in  material  between  the  hole  and  free  surface 
early  in  the  loading  history.  The  resulting  load  and  crack  growth  were 
suppressed  compared  to  that  observed  in  laboratory  specimens.  Plastic  hinge 
formation  is  not  observed  for  actual  laboratory  specimens  and  may  be 
attributed  to  the  two  dimension  idealization  of  a  three  dimensional  process  and 
mesh  size  used  in  that  region  of  the  model.) 

Material  Responses  -  Material  parameters  for  the  computational 
simulations  (Figure  2)  were  based  on  the  nonlinear  Cauchy  stress  versus 
logarithmic  ^-ain  curve  for  HY-100  steel.  The  response  of  HY-100  steel  is 
representative  of  extremely  ductile  steel  alloy  response.  The  HY-100  steel 
nonlinear  stress  versus  strain  parameters  (Figures  2c.  2d  and  2e)  used  in  the 
analysis  were  previously  obtained  with  demonstrated  accuracy  at  large  strains 
(Matic,  Kirby  and  Jolles,  [26]).  The  influence  of  material  strength  on  specimen 
strength  and  crack  growth  was  considered  by  introducing  lower  (Figure  2a 
and  2b)  and  higher  yield  stress  (Figure  2f  and  2g)  material  parameters.  The 
nonlinear  form  of  each  high  and  low  yield  stress  curve,  beyond  the  yield 
point,  followed  from  the  moderate  yield  HY-100  response.  The  HY-100  stress 
versus  strain  curves  exhibited  a  yield  stress  of  7.17  x  102  MN/m2  (1.04  x  10^ 
ib/in2).  The  high  and  low  yield  stress  curves  were  shifted  by  plus  and  minus 
1.37  x  102  MN/m2  (2.00  x  10^  lb/in2),  respectfully,  to  6.80  x  102  MN/m2  (0.84  x 
105  lb/in2)  and  8.54  x  102  MN/m2  (1.24  x  105  lb/in2). 

An  incremental  rate  independent  plasticity  theory  was  used  for  the 
material  constitutive  model  [27],  This  standard  model  for  plasticity  is 
summarized  here  for  completeness.  Total  strains  in  the  multiaxial  strain  state 
eij  were  obtained  by  the  integration  of  the  linearly  decomposed  rate  of 
deformation  tensor  Djj  This  integration  was  performed  under  the  assumption 
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Fig.  2  —  Cauchy  stress  versus  log  strain  curves  for  yield  stress  (  x  102  MN/m2)  —  fracture  toughness  (  x  102  MN-m/m3) 
combinations  (a)  5.79^.41,  (b)  5.79-13.23,  (c)  7.17-4.41,  (d)  7.17-8.83,  (e)  7.17-13.23,  (0  8.55-4.41,  (g)  8.55-13.23 
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Fig.  2  (Continued)  —  Cauchy  stress  versus  log  strain  curves  for  yield  stress  (  x  102  MN/m2)  —  fracture  toughness  (  x  102 
MN-m/m3)  combinations  (a)  5.79-1.41,  (b)  5.79-13.23,  (c)  7.17-4.41,  (d)  7.17-8.83,  (e)  7.17-13.23,  (0  8.55-4.41,  (g) 
8.55-13.23 
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Fig.  2  (Continued)  —  Cauchy  stress  versus  log  strain  curves  for  yield  stress  (  x  102  MN/m2)  —  fracture  toughness  (  x  102 
MN-m/m5)  combinations  (a)  5.79-4.41,  (b)  5.79-13.23,  (c)  7.17-4.41,  (d)  7.17-8.83,  (e)  7.17-13.23,  (0  8.55-4.41,  (g) 
8.55-13.23 
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Fig.  2  (Continued)  —  Cauchy  stress  versus  log  strain  curves  for  yield  stress  (  x  102  MN/mz)  —  fracture  toughness  (  x  102 
MN-m/m3)  combinations  (a)  5.79-4.41,  (b)  5.79-13.23,  (c)  7. 17-4.41,  (d)  7.17-8.83,  (e)  7.17-13.23,  (f)  8.55-4.41,  (g) 
8.55-13.23 


13 


200 


160 


n 

o 

r- 

X 

'b 


120 


co 

CO 

LU 

cc 

I— 

CO 

> 

X 

o 

< 

o 


80 

40 


Fig.  2  (Continued)  —  Cauchy  stress  versus  log  strain  curves  for  yield  stress  (  x  102  MN/m2)  —  fracture  toughness  (  x  102 
MN-tn/m3)  combinations  (a)  5.79-4.41,  (b)  5.79-13.23,  (c)  7.17-4.41,  (d)  7.17-8.83,  (e)  7.17-13.23,  (f)  8.55-4.41,  (g) 
8.55-13.23 
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Fig.  2  (Continued)  —  Cauchy  stress  versus  log  strain  curves  for  yield  stress  (  x  102  MN/m2)  —  fracture  toughness  (  x  102 
MN-m/m3)  combinations  (a)  5.79-4.41,  (b)  5.79-13.23,  (c)  7.17-4.41,  (d)  7.17-8.83,  (e)  7.17-13.23,  (0  8.55-4.41,  (g) 
8.55-13.23 
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Fig.  2  (Continued)  —  Cauchy  stress  versus  log  strain  curves  for  yield  stress  (  x  102  MN/m2)  —  fracture  toughness  (  x  102 
MN-m/nv*)  combinations  (a)  5.79-4.41,  (b)  5.79-13.23,  (c)  7.17-4.41,  (d)  7.17-8.83,  (e)  7.17-13.23,  (f)  8.55-4.41,  (g) 
8.55-13.23 
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that  the  elastic  strains  remain  infinitesimal,  as  are  the  cases  in  this 
investigation.  The  total  multiaxial  strain  state  eij,  expressed  in  terms  of  elastic 
and  plastic  components,  was 


+  £?. 
y 


(1) 


The  total  logarithmic  uniaxial  strain  £,  consistent  with  the  integration  of  the 
rate  of  deformation  tensor  for  a  multiaxial  strain  state,  was  decomposed  as 

£  =  £*  +  £?  (2) 

The  yield  function  f  takes  the  form 

f(Tij)  =  T(£p)  (3) 

where  tij  and  T  are  the  multiaxial  Kirchoff  (or  Treffetz)  and  uniaxial  stress 
states,  respectively.  The  associated  flow  rule  governed  plastic  strain 
increments  by  the  relation 


d£?=X— 
J  axij 


(4) 


In  the  case  of  purely  elastic  behavior  X  =  0 .  For  active  material  yielding 


X>0  (5) 


X  =  d£P— 1— 

3f 

Tija~ 

(6) 

Plastic  strain  increments  also  satisfied  a  dissipation  equivalence 
condition 


X  d£P  =  Xjj  d£jj 


and  a  consistency  condition 


(7) 


17 


(8) 


—  dtij  -  —  d£P  =  0 

3t;j  3eP 


The  von  Mises  yield  function 


f(^ij)  —  \  (SijSjj)2 


(9) 


was  employed.  The  deviatoric  stress  tensor  Sjj  was  defined  as 


Sij  —  Tij  *  ^'Xyc8ij 


(10) 


where  the  hydrostatic  component  of  stress  is  ^kk /3. 


The  Kirchoff  stress  and  logarithmic  strain  measures  are  employed 
because  of  advantages  gained  in  computational  implementation.  The  Kirchoff 
stress  tensor  *ij  is  approximately  equal  to  the  more  physically  motivated 
Cauchy  stress  tensor  <*ij  for  deformations  involving  only  small  changes  in 
volume.  This  condition  was  implicit  in  these  analyses.  The  uniaxial  Cauchy 
stress  -  log  strain  constitutive  response  of  the  material  was  formally  input,  in 
multilinear  form,  as  Cauchy  stress  and  logarithmic  strain  pairs  for  the  ABAQUS 
program. 


Local  Fracture  Toughness  Criterion  -  The  material  fracture  toughness, 
as  characterized  by  the  absorbed  strain  energy  density  at  fracture,  has  been 
determined  to  vary  by  at  least  fifty  percent  in  the  case  of  HY-100  steel  [26] 
subjected  to  different  local  stress  versus  strain  histories.  The  role  of  multiaxial 
loading  histories  affects  both  inelastic  deformation  and  fracture. 
Fractographic  dependence  for  HY-100  steel  has  been  examined  by  Harvey  and 
Jolles  [28]  who  found  significant  differences  in  relative  microvoid  and 
cleavage  proportions  on  cylinderical  tensile  specimen  fracture  surfaces. 
Correlation  between  these  microstructural  differences  and  the  absorbed  strain 
energy  density  at  the  point  of  initial  material  fracture  was  also  demonstrated. 
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In  order  to  examine  the  influence  of  a  local  fracture  toughness  value  on 
crack  growth,  the  path  dependent  aspects  of  fracture  toughness  were 
suppressed  by  considering  a  constant  critical  absorbed  strain  energy  density 
for  each  computational  simulation,  i.e. 

wc  =  constant.  (11) 

The  actual  path  dependent  of  the  toughness  was  used  to  select  the  range  of  wc 
values  for  the  simulations.  Three  different  critical  energy  density  values  were 
used.  The  low,  moderate  and  high  local  fracture  toughness  values  were  4.41, 

8.83  and  13.23  x  102  MN-m/m3  (64.0,  128.0  and  192.0  x  103  lb-in/in3).  The 
moderate  and  high  values  were  based  on  the  range  of  path  dependent  critical 
HY-100  values  obtained  from  the  solution  procedure  discussed  previously  [26]. 
The  importance  of  this  range  to  material  fracture  and  crack  initiation  of  HY- 

100  steel  compact  tension  specimens  has  also  been  previously  demonstrated 
DeGiorgi,  Kirby  and  Jolles  [18].  The  low  value  was  selected  to  simulate  what 
could  be  expected  if  composition  or  structure  changes  were  induced  in  the 
material  by  processes  such  as  welding. 

Computational  Model  -  The  ABAQUS  finite  element  code  Version  4.8  was 
used  for  the  computational  simulations  [29].  Finite  element  models  for  the 
simulations  are  shown  in  Figure  3.  Specimen  symmetry  required  that  only  one 
half  of  the  specimen  be  modeled.  The  initial  crack  faces  and  crack  propagation 
path  lie  on  the  symmetry  axis. 

Type  CPE8  and  CPE8H  plane  strain  elements  were  used  to  model  the  solid. 
Both  element  types  are  8-noded  quadrilateral  elements  with  quadratic 
displacement  interpolation.  The  CPE8H  elements  also  include  an  independently 
interpolated  linear  hydrostatic  stress.  The  CPE8H  elements  were  employed  in 

elements  near  crack  tip  to  sustain  the  anticipated  large  strain,  essentially 
constant  volume  deformation.  Interface  elements,  type  INTER3,  bonded  the 

material  elements  adjacent  to  the  Mode  I  crack  propagation  path  to  the  line  of 
symmetry.  The  interface  elements  were  the  means  by  which  nodal  force 
release  was  implemented. 


(a) 


(b) 

Fig.  3  -  Finite  element  models  for  ductile  tear  simulations  (a)  short  crack  length,  (b)  long  crack  length 
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All  analyses  were  performed  with  full  geometric  nonlinearity  to 

account  for  large  strains  and  large  rotations.  An  updated  Lagrangian 
formulation  was  used  for  incremental  solution  in  ABAQUS.  A  modified  Rik’s 
algorithm,  being  load  and  displacement  controlled,  was  used  to  enhance 
numerical  stability. 

Stress  boundary  conditions  were  applied  at  the  upper  left  portion  of  the 

model  to  approximate  the  stress  applied  at  a  loading  hole  in  an  actual  compact 
fracture  specimen.  More  detailed  consideration  of  the  loading  hole,  as 

considered  in  [18],  would  have  increased  the  model  size  but  not  provided  any 

additional  substantive  information  on  crack  growth,  crack  length  effects  or 

material  effects. 

Debonding  Algorithm  -  An  automated  approach  to  quasistatic  crack 
propagation  was  implemented  in  the  ABAQUS  Version  4.8  finite  element 

program.  The  approach  was  based  on  the  debonding  of  finite  elements, 
through  the  local  fracture  criterion  and  nodal  force  release,  to  create  new 
fracture  surface  during  the  load  application.  Nodal  release  was  implemented 
through  the  use  of  the  interface  elements  which  enforced  displacement 

continuity  across  a  specified  crack  propagation  path  ahead  of  the  crack  tip. 

The  crack  path  was  specified  in  the  model  by  a  list  of  nodes  ahead  of  the  crack 
tip.  Interface  elements  associated  with  the  crack  path  were  also  specified. 

A  subroutine,  written  by  the  user  and  executed  by  the  ABAQUS  program, 
was  used  to  determine  the  value  of  a  debond  parameter  c  at  each  node  of  the 
crack  path  ahead  of  the  current  crack  tip.  The  parameter  c,  defined  as  the 

difference  between  the  absorbed  strain  energy  density  w  and  associated 

critical  value  of  the  critical  strain  energy  density  wc,  were  associated  with 
each  node  in  the  crack  path  through  the  user  subroutine.  That  is, 

c  =  w  -  wc 

The  value  of  c  was  used  to  control  the  debonding  sequence  used  to 
predict  crack  propagation  along  the  specified  path  (Figure  4).  Fracture 

occurcd  when  the  debond  parameter  at  the  crack  tip  passed  from  negative  to 
positive  values,  reflecting  an  excess  of  absorbed  strain  energy  density  over 
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c  >  o  CONTINUE  DEBOND 


c  <  o  NO  DEBOND 


Fig.  4  —  Schematic  of  debonding  algorithm 
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the  material's  capacity  to  sustain  such  an  energy  density.  In  general,  a  load 

increment  did  not  produce  a  value  of  c  equal  to  zero.  Rather,  c  would  exceed 

zero  for  the  solution  increment  when  crack  growth  should  occur.  The 
debonding  algorithm  interpolated  over  the  normal  load  increment  to  generate 

an  intermediate  load  increment  leading  up  to  debonding  which  more  closely 
approximated  a  c  value  of  zero. 

Crack  propagation  occured  in  discrete  increments,  one  element  at  a 
time.  The  energy  density  parameters  were  calculated  during  each  ioad 

increment  from  the  material  constitutive  model  and  the  local  material  fracture 

toughness  criterion.  Required  information  was  obtained  from  the  finite 

element  database  by  the  user  subroutine. 

When  the  debonding  sequence  for  an  element  was  activated  the  load  was 
increased  by  a  small  fraction  of  the  current  load.  Therefore,  crack 
propagation  occured  under  essentially  constant  load.  The  two  nodes  of  the 

eight  noded  finite  elements  in  the  crack  path  immediately  ahead  of  the 

current  crack  tip,  i.e.  one  corner  and  one  midside  node,  were  released  by 
removing  displacement  constraints  and  reaction  forces.  The  displacement 
constraints  were  removed  immediately.  Nodal  forces  were  removed  in  equal 
increments.  Five  equal  increments  were  used  in  this  investigation. 
Equilibrium  was  reestablished  by  the  program  for  the  new  crack  length.  The 
fracture  criterion  was  examined  to  see  if  crack  growth  should  continue  at  the 

current  load.  If  the  fracture  criterion  was  satisfied,  crack  growth  continues  at 
the  same  step  in  the  load  history.  If  the  criterion  was  not  satisfied,  loading 
proceeds  to  the  next  increment.  Stable  crack  growth  was  viewed  as  a  process  of 
continuous  loading  accompanied  by  intermittent  debonding.  Unstable  crack 

growth,  if  it  occurred,  would  be  viewed  as  continuous  debonding  in  the 
absence  of  additional  active  loading. 

The  simulations  conducted  for  both  long  and  short  crack  length 
geometries  are  summarized  in  Table  2.  The  extreme  values  of  all  strength  and 

local  fracture  toughness  combinations  were  examined,  as  were  all  local 
fracture  toughness  values  for  the  HY-100  constitutive  response.  These 

simulations  were  identified  by  an  abbreviation  of  the  form  GMMFFF  and 
associated  symbols  in  the  graphical  results.  G  denotes  the  short  (S)  or  long  (L) 
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Yield  Stress,  Local  Fracture  Toughness,  x  10^  MN-m/m^ 

x  lO^  MN/m^  Low  Moderate  High 


4.41  8.83  13.23 


Table  2  Material  Parameter  Variations  Analyzed  for  Short 
and  Long  Crack  Length  Geometries 
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crack  geometry  and  is  represented  by  open  and  closed  symbols,  respectively, 
on  all  graphs.  MM  is  the  material  yield  strength.  Low  (-2),  moderate  (HY),  and 
high  (+2)  yield  stresses  are  represented  by  circles,  triangles  and  squares.  FFF 
is  the  local  fracture  toughness.  Low  (064),  moderate  (128)  or  high  (192)  critical 
absorbed  energy  density  parameters  are  indicated  by  dotted,  dashed  and  solid 
lines,  respectively. 

The  debond  criterion  was  explicitly  based  on  nodal  averaged  values  of 
energy  density  calculated  from  adjacent  element  nodal  averages.  As  is  common 
in  high  gradient  regions  of  a  finite  element  analysis,  the  nodal  averaged  value 
may  not  be  an  accurate  indication  of  material  response  immediately  ahead  of 
the  crack.  Specifically,  the  nodal  average  will  be  greater  than  the  value 
behind  the  crack  and  less  than  the  value  ahead  of  the  crack.  To  achieve  a  more 
consistent  interpretation  of  the  nodal  averaged  value  for  debonding,  a 
proportionality  factor  was  used  to  relate  the  nodal  averaged  energy  density 
values  to  the  energy  density  values  of  the  material  in  the  element  ahead  of  the 
crack  tip.  This  factor,  based  on  data  in  previous  studies  [18],  reduced  the  low, 
moderate  and  high  critical  energy  density  values  to  2.41,  4.83  and  7.23  x  102 
MN-m/m^  (35.0,  70.0  and  105.0  x  10^  lb-in/in^),  respectively,  for  purposes  of 
numerical  implementation  in  the  debonding  algorithm.  The  physical  response 
of  the  material  and  the  fracture  toughness  value  are  therefore  both 
considered  by  the  proportionality  factor  in  a  consistent  manner. 


V  Results 

The  results  of  the  crack  growth  simulations  are  presented  here  in  terms 
of  local  material  responses  and  global  specimen  responses. 

The  strain  energy  density  ratio  (w/wc)  at  the  current  crack  tip  is 
plotted  against  crack  mouth  opening  displacement  (CMOD).  For  the  three 
materials  involving  the  moderate  yield  strength  parameters  associated  with 
HY-100  steel  the  short  crack  length  results  are  shown  in  Figure  5  and  the  long 
crack  length  results  are  shown  in  Figure  6.  The  four  extreme  yield  strength 
and  local  fracture  toughness  parameter  combinations  for  the  short  crack 
length  and  the  long  crack  length  are  shown  in  Figure  7  and  Figure  8. 
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Fig.  5  Strain  energy  density  ratio  (w/wc)  at  the  crack  tip  versus  crack  mouth  opening  displacement  (CMOD)  for  HY-100 
short  crack  length,  7.17  x  102  MN/m2  ksi  yield  stress,  and  (a)  4.41,  (b)  8.83  and  (c)  13.23  x  102  MN-m/m’  local 
fracture  toughness  values 
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Fig.  6  —  Strain  energy  density  ratio  (w/wc)  at  the  crack  lip  versus  crack  mouth  opening  displacement  (CMOD)  for  HY-100 
long  crack  length.  7.17  x  102  MN/m2  yield  stress,  and  (a)  4.41,  (b)  8.83  and  (c)  13.23  x  I02  MN-m/m5  local  fracture 
toughness  values 
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Fig.  7  —  Strain  energy  density  ratio  (w/wc)  at  the  crack  tip  versus  crack  mouth  opening  displacement  (C MOD)  for  short 
crack  length  geometry  and  extreme  combinations  of  strength  and  toughness.  Yield  stress  of  5.79  x  102  MN/m2  and  (a) 
4.41  and  (b)  13.23  x  102  MN-m/m3  local  fracture  toughness  values.  Yield  stress  of  8.55  x  102  MN/m2  and  (c)  4.41  and 
(d)  13.23  x  102  MN-m/m3  local  fracture  toughness  values. 
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Fig.  7  (Continued)  —  Strain  energy  density  ratio  (w/wc)  at  the  crack  tip  versus  crack  mouth  opening  displacement  (CMOD) 
for  short  crack  length  geometry  and  extreme  combinations  of  strength  and  toughness.  Yield  stress  of  5.79  x  102  MN/m2 
and  (a)  4.41  and  (b)  13.23  x  102  MN-mim3  local  fracture  toughness  values.  Yield  stress  of  8.55  x  102  MN/m2  and  (c) 
4.41  and  (d)  13.23  x  102  MN-m/m3  local  fracture  toughness  values. 
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Fig.  8  —  Strain  energy  density  ratio  (w/wc)  at  the  crack  tip  versus  crack  mouth  opening  displacement  (CMOD)  for  long 
crack  length  geometry  and  extreme  combinations  of  strength  and  toughness.  Yield  stress  of  5.79  x  102  MN/m1  and  (a) 
4.41  and  (b)  13.23  x  102  MN-m/m’  local  fracture  toughness  values.  Yield  stress  of  8.55  x  102  MN/m2  and  (c)  4.41  and 
(d)  13.23  x  102  MN-m/m3  local  fracture  toughness  values. 
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Fig.  8  —  (Continued)  Strain  energy  density  ratio  (w/wc)  at  the  crack  tip  versus  crack  mouth  opening  displacement  (CMOD) 
for  long  crack  length  geometry  and  extreme  combinations  of  strength  and  toughness.  Yield  stress  of  5.79  x  102  MN/m2 
and  (a)  4.41  and  (b)  13.23  x  102  MN-m/m3  local  fracture  toughness  values.  Yield  stress  of  8.55  x  102  MN/m2  and  (c) 
4.41  and  (d)  13.23  x  102  MN-m/m3  local  fracture  toughness  values. 
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The  strain  energy  density  ratio  normalizes  the  material  damage  with 

respect  to  the  material  resistance  to  fracture  at  any  point  in  the  specimen  over 
the  loading  history.  The  rate  at  which  a  point  in  the  specimen  approaches 
fracture  is  reflected  in  the  slope  of  the  curve.  A  strain  energy  density  ratio  of 
one  corresponds  to  material  fracture  and  initiation  of  the  debonding  process. 

The  current  crack  tip  position  is  redefined  by  the  debonding  process 

when  crack  growth  takes  place.  The  energy  density  ratio  of  the  material  at  the 
new  crack  tip  again  represents  the  local  resistance  to  fracture  .The  new  crack 
tip  must  pass  through  this  material  if  crack  growth  is  to  be  sustained.  For  the 
computational  simulations  of  this  study,  the  decrease  in  strain  energy  density 
ratio  observed  after  each  successive  crack  growth  increment  of  one  element 
results  in  a  new  stable  crack  length.  For  each  successive  crack  growth 
increment,  however,  the  energy  density  ratio  of  the  material  at  the  new  crack 
tip  is  greater  than  the  energy  density  ratio  at  the  time  the  prior  crack  tip  was 
defined  by  debonding.  This  trend,  as  defined  by  the  energy  density  ratios  at 
successive  crack  tip  definitions,  is  an  indication  of  the  rate  at  which  stable 
crack  growth  leads  to  crack  instability.  For  those  material  parameters 
considered  for  which  multiple  debonding  increments  occurred,  it  is  apparent 
that  more  than  ten  total  crack  growth  increments  would  be  required,  at  the 

present  rate  of  post-increment  energy  ratio  increase,  to  produce  repetitive 
debonding  and  crack  instability.  More  than  ten  increments  of  crack  growth 
would  place  the  crack  tip  well  into  the  center  portion  of  the  specimen  width. 
The  range  of  parameter  values  would  generally  suggest  this  for  all  cases 
considered. 

Physically,  HY-100  steel  in  compact  fracture  specimen  geometries  is 
observed  to  sustain  extensive  stable  crack  growth  at  the  leading  edge  of  the 
crack,  in  the  center  of  the  specimen,  extending  over  more  than  one  half  the 
specimen  width.  The  crack  growth  at  the  free  surface  lags  far  behind, 
resulting  in  the  classical  curved  crack  front  which  tunnels  into  the  specimen. 

(Plastic  hinge  formation  occurs  as  the  crack  approaches  the  free  surface  of 
the  specimen  ahead  of  the  crack.  Final  specimen  failure  follows.)  The  physical 
response  is  an  inherently  three  dimensional  process,  while  the  simulations 
were  two  dimensional.  The  nature  of  these  results  are,  however,  consistent 
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with  the  general  trends  of  HY-100  crack  growth  preceding  the  gross 
geometric  deformations  of  plastic  hinge  formation. 

The  applied  force  is  plotted  against  crack  mouth  opening  displacement 
(CMOD),  taken  at  the  comer  defined  by  the  intersection  of  the  crack  face  and 
free  surface  behind  the  crack  on  the  specimen.  For  the  simulations  involving 
HY-100  yield  strength  the  short  crack  length  results  are  shown  in  Figure  9 
and  the  long  crack  length  results  are  shown  in  Figure  10.  The  four  extreme 
yield  strength  and  local  fracture  toughness  combinations  for  the  short  crack 
length  and  the  long  crack  length  are  shown  in  Figure  11  and  Figure  12.  The 

overall  nonlinearity  of  the  response  is  sustained  though  increasing  values  of 

CMOD.  The  intermittent  unloading  on  the  curves,  interrupting  the  otherwise 
smooth  character  of  the  response,  are  generated  as  debonding  extends  the 
crack.  The  unloading  occurs  during  the  debond  as  equilibrium  is  reestablished 
by  the  modified  Rik's  algorithm  for  a  longer  crack  and  corresponding  elastic- 

plastic  crack  tip  fields.  Active  loading  of  the  specimen  resumes  after 
equilibrium  is  obtained.  Relatively  small  unloading  magnitudes  encountered 
during  stable  debonding  suggests  that  the  discrete  debonding  is  approximating 
continuous  stable  crack  growth  in  a  satisfactory  manner. 

It  may  be  noted  that  the  influence  of  greater  material  yield  stress  is 
directly  translated  into  higher  specimen  strength.  At  the  same  time,  crack 
growth  rates  are  seen  to  accelerate  as  the  yield  stress  increases  for  constant 

material  fracture  toughness  and  specimen  geometry.  The  elastic  portion  of  the 
material  response,  locally,  and  its  interaction  with  the  aggregate  recoverable 
response  of  the  specimen,  globally,  combine  to  produce  conditions  favorable  to 

energy  dissipation  at  the  crack  tip.  This  apparently  accelerates  the  local 

fracture  process  and,  therefore,  crack  growth. 

It  may  also  be  noted  that  the  influence  of  material  fracture  toughness 
on  specimen  strength,  as  seen  for  the  case  of  HY-100  material  yield  strength, 

was  most  significant  when  comparing  the  low  to  moderate  fracture  toughness 
simulations.  The  relative  ability  of  the  specimen  to  sustain  load  is  increased 
only  slightly  when  the  material  fracture  toughness  was  increased  from  the 
moderate  to  high  value.  A  comparison  of  the  uniaxial  stress  and  strain  ranges 
associated  with  the  local  fracture  toughness  values  (i.e.  Figure  2c)  shows  that 
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0.5*CMOD  (X  10"2  CM) 

Fig.  9  —  Applied  load  (P)  versus  crack  mouth  opening  displacement  (CMOD)  for  HY-100.  Short  crack  length,  7.17  x  102 
MN/m2  yield  stress  and  4.41,  8.83  and  13.23  x  102  MN-m/m3  local  fracture  toughness  values 
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Fig  10  —  Applied  load  (P)  versus  crack  mouth  opening  displacement  (CMOD)  for  HY-100.  Long  crack  length,  7.17  x 
10*  MN/nr  yield  stress  and  4.41,  8.83  and  13.23  x  102  MN-m/m3  local  fracture  toughness  values. 
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Fig.  12  —  Applied  load  (P)  versus  crack  mouth  opening  displacement  (CMOD)  for  extreme  combinations  of  strength  and 
toughness.  Long  crack  length,  5.79  and  8.55  x  102  MN/m2  yield  stress,  and  4.41  and  13.23  x  102  MN-m/m}  local 
fracture  toughness  values. 
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the  low  fracture  toughness  value  includes  only  the  pronounced  strain 
hardening  portion  of  the  HY-100  nonlinear  constitutive  response.  The 
moderate  and  high  toughness  values  (i.e.  Figures  2d  and  2e)  include  both  the 
pronounced  and  weak  strain  hardening  portions  of  the  HY-100  nonlinear 
constitutive  response.  The  moderate  and  high  toughness  materials  are 
apparently  able  to  sustain  the  additional  dissipated  energy  without  a 
significant  increase  in  elastic  energy  associated  with  this  dissipation.  The 

local  elastic  energy  is  also  distributed  in  a  more  uniform  manner  over  regions 
of  large  effective  plastic  strain.  This  apparently  allows  for  more  distributed 

local  loading  and  unloading  associated  with  crack  growth. 

Each  simulation  was  planned  to  run  for  sixty  load  increments.  Since  the 
debonding  sequence  was  not  specified  a  priori  and  a  range  of  material 
parameters  were  used,  the  amount  of  crack  growth  varied  considerably.  Crack 
growth  for  long  and  short  crack  length  simulations  are  given  in  Tables  3  and 

4,  respectively.  In  several  instances,  combinations  of  material  parameters, 
geometry,  load  and  mesh  were  unable  to  maintain  convergence  with  further 
loading  or  debonding  at  some  point  in  the  analysis.  Material  stress  and  strain 

limited  convergence  during  crack  growth  was  noted  for  some  of  the 
simulations.  These  occurrences,  their  cause  and  solution  are  discussed  here. 

For  the  case  of  high  yield  stress,  elastic  stresses  increase  with  yield 
stress  for  a  given  value  of  equivalent  plastic  strain.  The  crack  growth 
increments  produce  simultaneous  global  unloading  and  increased  plastic 
deformation  at  the  new  crack  tip.  The  global  elastic  unloading  which 
developed  during  an  increment  of  crack  growth  produced  large  displacements. 
This  reflected  the  decrease  in  specimen  compliance  associated  with  an 
increasing  crack  length  and  the  shift  of  the  crack  tip  plastic  zone. 

High  toughness  critical  absorbed  strain  energy  density  values  were  also 
associated  with  increasing  computational  demands  on  the  simulations.  The 
relatively  low  slope  of  the  material  responses  produces  large  plastic  strains. 
The  coupled  effect  of  high  toughness  and  low  yield  stress  material  response  is 
capable  of  producing  the  largest  strains.  These  strains  became  difficult  to 
sustain  in  conjunction  with  the  debonding  for  the  longer  crack  simulations. 
Additional  refinement  in  mesh  size  would  facilitate  additional  debonding. 
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Yield  Stress, 

Local 

Fracture  Toughness,  x  102 

MN-m/m^ 

x  102  MN/tn2 

Low 

Moderate 

High 

4.41 

8.83 

13.23 

High 

8.55 

4 

“ 

3 

Moderate 

7.17 

4 

3 

2 

Low 

5.79 

3 

2 

Elements 

a/ao 

1 

1.17 

2 

1.33 

3 

1.50 

4 

1.67 

Table  3  Number  of  elements  sustaining  crack  growth 
for  short  crack  length  geometry 
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Yield  Stress, 

Local  Fracture  Toughness,  x  102 

MN-m/m^ 

x  102  MN/m2 

Low 

Moderate 

High 

4.41 

8.83 

13.23 

High 

8.55 

1 

- 

2 

Moderate 

7.17 

4 

2 

2 

Low 

5.79 

3 

1 

Elements 

a/ao 

1 

1.09 

2 

1.18 

3 

1.27 

4 

1.36 

Table  4  Number  of  elements  sustaining  crack  growth 
for  long  crack  length  geometry 
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From  each  of  these  considerations,  the  following  conclusions  can  be 
drawn.  First,  for  simulations  of  relatively  short  length  cracks  sustaining 
significant  crack  growth,  the  global  load  versus  displacement  response  was 

piecewise  continuous  with  relatively  minor  discontinuities  reflecting  the 
discrete  debonding  process.  Additional  mesh  refinement  would  further  reduce 
these  discontinuities  by  providing  finer  crack  increments.  This  modification 
is  not  necessary,  however,  to  assess  the  load  versus  displacement  response  of 
the  specimen.  Second,  large  strains  associated  with  low  yield  stress  and  high 
local  fracture  toughness  reached  limits  governed  by  the  mesh  size  ahead  of 
the  crack  tip  utilized  in  this  investigation.  Additional  mesh  refinement  would 
be  necessary  to  sustain  crack  growth  preceded  by  effective  log  plastic  strains 

over  1.0.  Third,  the  highest  yield  stress  materials  considered  generated 
significant  elastic  displacements  as  a  result  of  unloading  which  occured 
during  debonding  increments.  In  this  case,  additional  mesh  refinement  would 
reduce  the  magnitude  of  unloading  sustained  by  each  debond  increment. 

It  should  be  noted  that  many  materials  have  lower  yield  stresses,  sustain 
smaller  strains  and  respond  with  a  lower  local  fracture  toughness  than  HY-100 
or  the  parametric  variations  considered  here.  In  view  of  the  relative  size  of 
the  computational  resources  allocated  to  any  individual  simulation,  the  ability 

to  conduct  analyses  based  on  approaches  such  as  those  presented  here  are 
currently  feasible  and  will  continue  to  grow  as  improvements  in  computer 

technology  continue  to  become  available. 


VI  Summary 

Parametric  computational  studies  of  stable  crack  growth  in  a  modified 
compact  tension  specimen  geometry  were  conducted.  Computational  modeling 
issues  of  crack  growth  prediction  using  a  local  fracture  criterion  were 
addressed  using  finite  element  analysis  incorporating  a  debonding  algorithm. 
The  debonding  algorithm  utilized  a  user  subroutine  to  assess  local  fracture  at 
the  crack  tip  and  a  nodal  release  scheme  to  produce  crack  growth.  Crack 
growth  was  treated  as  a  dependent  variable  of  the  analysis. 
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HY-100  steel  material  parameters  were  used  as  a  basis  for  nonlinear 
constitutive  response  and  local  fracture  toughness.  Reduced  and  enhanced 
material  yield  strength  and  fracture  toughness  were  examined  to  assess  their 
relative  effects  on  local  fracture  at  the  crack  tip,  the  rate  at  which  the  crack 
stability  is  approached  and  the  global  nonlinearity  produced  by  the 
combination  of  inelastic  material  response  and  crack  growth.  A  comparison  of 
constant  yield  strength  and  local  fracture  toughness  simulations  suggest  that 
while  stress  and  strain  levels  which  a  material  can  sustain  are  important,  the 
manner  in  which  these  stress  and  strain  levels  are  obtained  through 
constitutive  nonlinearity  is  equally  important  to  specimen  strength  and  crack 
growth  response. 

The  combined  influence  of  the  material  responses  and  the  finite 
element  model  on  solution  convergence  for  this  type  of  analysis  was  discussed. 
The  relative  influence  of  high  elastic  stresses  and  large  plastic  strain  on  both 
local  and  global  unloading  encountered  during  nodal  release  were  noted. 

The  results  of  this  study  suggest  that  the  capability  to  model  and 
simulate  stable  crack  growth  in  two  dimensions  using  local  fracture  criteria  is 
computationally  accessible.  The  simulations  conducted  in  this  investigation, 
representing  strong,  ductile  and  fracture  resistant  alloys,  encompass  the 
computational  modeling  requirements  of  many  materials  currently  in  use.  The 
capability  to  extend  local  fracture  criteria  and  crack  growth  simulation  into 
three  dimensions  will  be  realized  as  computational  technology  continues  to 
evolve. 
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